########################################################################
############################ R CODE TO RUN  ############################
####################### ECOLOGICAL NICHE MODELS ########################
############################################################# March/2024

## Dr.Neander M. Heming - PPGECB - UESC
## www.researchgate.net/profile/Neander_Heming
## neanderh@yahoo.com.br
##===========================================================##
## 
## 

{
  library(rJava)
  library(dismo)
  maxent()  
  library(readxl)
  library(ENMwizard)
}

# list.files("data/shapes/")

# occ new
occ.new <- shapefile("data/shapes/Pontos_Novos_Xanthosternos/Pontos_Sxanthosternos.shp")
hydro <- shapefile("data/shapes/Hidrografia/Hidrografia_area_Waldney.shp")
hydro2 <- shapefile("data/shapes/Hidrografia/HIDRO.shp")
AF <- shapefile("data/shapes/mata_atlantica/dominio_merge2.shp")

# Load occ data
spp.occ.list <- list()
spp.occ.list$Sapajus_xanthosternos <- shapefile("data/shapes/Sapajus_xanthosternos.shp")


#####------- 1. Prepare environmental data
# 1. Prepare enviromental layers
##### 1.1 create occ polygon to crop rasters prior to modelling
occ.polys <- set_calibarea_b(spp.occ.list, plot=F)

# 1.2 creating buffer
occ.b <- buffer_b(occ.polys, width = 1.5, plot=F)

# 1.3. Cut enviromental layers with M and save in hardrive.
# biovars.all <- sort(unique(c(c(1,3,4,5,7,12,15,19), # prim
#              c(1,3,5,12,15,16,18), # Bradypus
#              c(2,3,5,7,14,16,19)))) # Chaetomys
# predictors <- brick("data/rasters/envData.Sapajus_xanthosternos.gri") # PPI_AF spp.occ.list
predictors <- stack(brick("data/rasters/wc2.1_bio/wc2.1_30s_bio_AF.grd"),
                    raster("data/rasters/DrySeason100AF/DrySeason100_current_AF.grd")) # PPI_AF spp.occ.list
names(predictors)
predictors <- predictors[[c(paste0("bio", c(1,3,4,5,7,12,15,17,19)),"DrySeason100")]]
occ.b.env <- cut_calibarea_b(occ.b, predictors)
plot(occ.b.env$Sapajus_xanthosternos)

# 1.4 remove correlated variables from our variable set
vars <- select_vars_b(occ.b.env, cutoff=.75, sample.size=10000, names.only =T)
occ.b.env <- select_vars_b(occ.b.env, vars, cutoff=.75, names=F)

#####------- 2. Filtering original dataset # clean occ data
# {
#   occ.locs2 <- thin_b(spp.occ.list, thin.par = 3, reps = 3)
#   occ.locs2 <- load_thin_occ(occ.locs2)
#   # i <- names(occ.locs)[1]
#   for(i in names(occ.locs)){
#     occ.locs2[[i]] <- occ.locs[[i]]
#   }
#   occ.locs <- occ.locs2
#   rm(occ.locs2)
# }
occ.locs <- thin_b(spp.occ.list, thin.par = 10, reps = 3)
occ.locs <- load_thin_occ(occ.locs)
cbind(Original = sapply(spp.occ.list, nrow), Filtrado=sapply(occ.locs, nrow))



#####

# ENMeval_res.lst <- readRDS(file = "3_out.MaxEnt/ENMeval_res.lst_2021-04-14")
method <- ifelse(sapply(occ.locs, nrow)<=15, "jackknife", "block")
ENMeval_res.lst <- ENMevaluate_b(occ.locs, occ.b.env,
                                 RMvalues = c(.5, 1, 1.5, 2.5, 3.5, 4.5, 5.5),
                                 ### full list of models
                                 fc = c("L", "P", "Q", "H", "T", 
                                        "LP", "LQ", "LH", "LT", "PQ", "PH", "PT", "QH", "QT", "HT"),
                                        # "LPQ", "LPH", "LPT", "LQH", "LQT", "LHT",
                                        # "PQH", "PQT", "PHT",
                                        # "LPQH", "LPQT", "LPHT", "LQHT", "PQHT", 
                                        # "LPQHT"),
                                 method = method,
                                 parallel = T, numCores = 5, algorithm = 'maxent.jar')

# 4.3 Run top corresponding models
mxnt.mdls.preds.cf <- ENMwizard::calib_mdl_b(ENMeval.o.l = ENMeval_res.lst, 
                                             a.calib.l = occ.b.env, 
                                             occ.l = occ.locs,
                                             mSel = "EBPM", dAICc = 2, numCores = 6)

saveRDS(ENMeval_res.lst, file = paste0("3_out.MaxEnt/ENMeval_res.lst_thin10", Sys.Date(), ".rds"))


#### 4.8 projection on climate scenarios
{
  ## climate scenarios
  proj.scn <- list(ncurrent = predictors)
  ## polygon of selected area
  poly.projection <- set_projarea_b(occ.polys, mult = .55, buffer=FALSE)#
  plot(poly.projection$Sapajus_xanthosternos)
  plot(spp.occ.list$Sapajus_xanthosternos, add=T)
  plot(occ.new, col="red", pch=19, add=T)
  plot(hydro, add=T)
  # plot(hydro2, add=T)
  
  ## cutting projection area
  proj.area <- cut_projarea_mscn_b(poly.projection, proj.scn)
}

## projecting
mxnt.mdls.preds.cf <- proj_mdl_b(mxnt.mdls.preds.cf, 
                                 a.proj.l = proj.area,
                                 numCores = 4, parallelTunning=T)

## apply threshold
mods.thrshld.lst <- thrshld_b(mxnt.mdls.preds.cf, thrshld.i = 4)

## plots
occ.locs.spdf <- occ.locs$Sapajus_xanthosternos
coordinates(occ.locs.spdf) <- ~Longitude+Latitude

plot(mods.thrshld.lst$Sapajus_xanthosternos$ncurrent$binary$mtp)
plot(spp.occ.list$Sapajus_xanthosternos, add=T)
plot(occ.new, col="red", pch=19, add=T)
plot(hydro, col="blue", add=T)
plot(AF, add=T)

pdf("out/Sapajus_xanthosternos_suitability.pdf")
plot(mods.thrshld.lst$Sapajus_xanthosternos$ncurrent$continuous$mtp)
plot(spp.occ.list$Sapajus_xanthosternos, add=T)
plot(occ.new, col="red", pch=19, add=T)
plot(hydro, col="blue", add=T)
plot(AF, border="forestgreen", add=T)
dev.off()

extract(mods.thrshld.lst$Sapajus_xanthosternos$ncurrent$continuous$mtp, occ.new)

dir.create("out")
raster::writeRaster(mods.thrshld.lst$Sapajus_xanthosternos$ncurrent$continuous$mtp, 
                    'out/Sapajus_xanthosternos_suitability',
                    format="GTiff")


